function [sigmult,penalt] = mhcov(para,data)
% [sigmult,PENALT] = MHCOV(PMOD,DATA)
%
% computes square root of covariance matrix to be used for proposing step 
% in Metropolis-Hastings algorithm.
%
% INPUT
%	PMOD: posterior mode
%   DATA: data
%
% OUTPUT
%	SIGMULT : square root of inverse Hessian
%   PENALT  : penalty at the posterior mode (optional)
% 
% generalized inverse procedure is taken from PINV.M

% Sungbae An
% created: 06/29/2004
% updated: 07/21/2005 

global prior_

%--------------------------------------------------------------------
% evaluate Hessian
%--------------------------------------------------------------------
hssn = nhess('rsfnmh',para,data);


% take negative of Hessian
X=-hssn;

%--------------------------------------------------------------------
% Variance: Pseudo-inverse of the negative of Hessian
%--------------------------------------------------------------------

% singular value decomposition
[U,S,V] = svd(X,0);
s = diag(S);

% rank of Hessian (theoretical)
r1 = sum(prior_.maskinv);

% evaluate rank
tol = size(X,1)*max(s)*eps;
r2 = sum(s > tol);

r = r1;		% use theoretical rank

% inverse of negative Hessian
S(:) = 0;
if (r == 0)
	error('every singular value is zero.');
else
    S(1:r,1:r) = diag(ones(r,1)./s(1:r));
	sigmult    = U*sqrt(S);							% sigmult
	logdetS    = sum(log(diag(S(1:r,1:r))));
	penalt     = (r/2)*log(2*pi) + 0.5*logdetS;		% posterior mode penalty
end

